#############################################################################################
# NKPH models
# define model-specific mappings and derivative functions
#############################################################################################


using DataFrames


#############################################################################################
#############################################################################################
# baseline model
# no demand factor reaction to macro variables
# demand factor reaction to short rate
# Treasury and Risky demand factors identical parameterization
# alpha and alpha_til identical parameterization
# siginv and rho_pi fixed

function map_params_semiC!(nkphsolvr::NKPHSolver)::Nothing
    # update all model objects
    # pull out objects
    M, Upsilon, J, N_params = nkphsolvr.M, nkphsolvr.Upsilon, nkphsolvr.J, nkphsolvr.N_params
    sdyn = nkphsolvr.sdyn
    eta0, eta1 = nkphsolvr.eta0, nkphsolvr.eta1
    sigma, a = nkphsolvr.sigma, nkphsolvr.a
    alpha0, alpha0_til = nkphsolvr.alpha0, nkphsolvr.alpha0_til
    AR_hat = nkphsolvr.AR_hat
    aux_macro_params = nkphsolvr.aux_macro_params

    # pull out individual params from x_dict
    x_dict = nkphsolvr.x_dict

    # update model objects (inplace)
    # direct updates
    alpha0[1] = x_dict[:a0]
    alpha0_til[1] = x_dict[:a0]
    eta0[1] = x_dict[:n0]

    sigma[1, 1] = x_dict[:s_i]
    sigma[2, 2] = x_dict[:s_d]
    sigma[3, 3] = x_dict[:s_zpi]
    sigma[4, 4] = x_dict[:s_zx]
    sigma[5, 5] = sigma[6, 6] = sigma[7, 7] = x_dict[:s_b]

    # dynamics
    # Upsilon matrix
    # i_t
    Upsilon[1, 1] = x_dict[:k_i]
    Upsilon[1, J+1] = -x_dict[:k_i] * x_dict[:p_ipi]
    Upsilon[1, J+2] = -x_dict[:k_i] * x_dict[:p_ix]
    # d_t
    Upsilon[2, 2] = x_dict[:k_d]
    Upsilon[2, J+1] = -x_dict[:k_d] * x_dict[:p_dpi]
    Upsilon[2, J+2] = -x_dict[:k_d] * x_dict[:p_dx]
    # z_pit
    Upsilon[3, 3] = x_dict[:k_zpi]
    # z_xt
    Upsilon[4, 4] = x_dict[:k_zx]
    # b_st/b_lt/btil_t
    Upsilon[5, 5] = Upsilon[6, 6] = Upsilon[7, 7] = x_dict[:k_b]
    #Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:k_b] * x_dict[:p_bi]
    Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:p_bi]
    # pi_t
    Upsilon[J+1, 3] = 1.
    Upsilon[J+1, J+1] = -aux_macro_params[:rho_pi]
    Upsilon[J+1, J+2] = x_dict[:delta_x]
    # x_t
    Upsilon[J+2, 1:J] = -aux_macro_params[:siginv] .* AR_hat
    Upsilon[J+2, 4] += aux_macro_params[:siginv]

    return nothing
end

# derivs wrt params
function deriv_param_semiC!(nkphsolvr::NKPHSolver, param_idx::Int)::Nothing
    # compute model object derivs wrt param_{param_idx}
    J, N_params = nkphsolvr.J, nkphsolvr.N_params
    # pull out deriv objects (views)
    dsigma = @view nkphsolvr.dsigma[:, :, param_idx]
    dUpsilon = @view nkphsolvr.dUpsilon[:, :, param_idx]
    dalpha0 = @view nkphsolvr.dalpha0[param_idx]
    dalpha0_til = @view nkphsolvr.dalpha0_til[param_idx]
    deta0 = @view nkphsolvr.deta0[param_idx]

    # deriv mapping depends on param_idx
    pmap = nkphsolvr.param_model_map[param_idx]
    pmap.vals[:] .= false

    # get param symbol
    param_sym = nkphsolvr.x_dict.keys[param_idx]


    ##### sigma params
    if param_sym == :s_i
        pmap[:sigma] = true
        dsigma[:, :] .= 0
        dsigma[1, 1] = 1.
    elseif param_sym == :s_d
        pmap[:sigma] = true
        dsigma[:, :] .= 0
        dsigma[2, 2] = 1.
    elseif param_sym == :s_zpi
        pmap[:sigma] = true
        dsigma[:, :] .= 0
        dsigma[3, 3] = 1.
    elseif param_sym == :s_zx
        pmap[:sigma] = true
        dsigma[:, :] .= 0
        dsigma[4, 4] = 1.
    elseif param_sym == :s_b
        pmap[:sigma] = true
        dsigma[5, 5] = dsigma[6, 6] = dsigma[7, 7] = 1.
    ##### Upsilon params: mean reversion
    elseif param_sym == :k_i
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        phi_pi, phi_x = nkphsolvr.x_dict[:p_ipi], nkphsolvr.x_dict[:p_ix]
        dUpsilon[1, 1] = 1.
        dUpsilon[1, J+1] = -phi_pi
        dUpsilon[1, J+2] = -phi_x
    elseif param_sym == :k_d
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        phi_pi, phi_x = nkphsolvr.x_dict[:p_dpi], nkphsolvr.x_dict[:p_dx]
        dUpsilon[2, 2] = 1.
        dUpsilon[2, J+1] = -phi_pi
        dUpsilon[2, J+2] = -phi_x
    elseif param_sym == :k_zpi
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        dUpsilon[3, 3] = 1.
    elseif param_sym == :k_zx
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        dUpsilon[4, 4] = 1.
    elseif param_sym == :k_b
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        #phi_i = nkphsolvr.x_dict[:p_bi]
        dUpsilon[5, 5] = dUpsilon[6, 6] = dUpsilon[7, 7] = 1.
        #dUpsilon[5, 1] = dUpsilon[6, 1] = dUpsilon[7, 1] = -phi_i
    ##### Upsilon params: macro reaction
    elseif param_sym == :p_ipi
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        kappa = nkphsolvr.x_dict[:k_i]
        dUpsilon[1, J+1] = -kappa
    elseif param_sym == :p_ix
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        kappa = nkphsolvr.x_dict[:k_i]
        dUpsilon[1, J+2] = -kappa
    elseif param_sym == :p_dpi
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        kappa = nkphsolvr.x_dict[:k_d]
        dUpsilon[2, J+1] = -kappa
    elseif param_sym == :p_dx
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        kappa = nkphsolvr.x_dict[:k_d]
        dUpsilon[2, J+2] = -kappa
    elseif param_sym == :p_bi
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        #kappa = nkphsolvr.x_dict[:k_b]
        #dUpsilon[5, 1] = dUpsilon[6, 1] = dUpsilon[7, 1] = -kappa
        dUpsilon[5, 1] = dUpsilon[6, 1] = dUpsilon[7, 1] = -1.0
    ##### Upsilon params: IS/PC macro terms (delta_x)
    elseif param_sym == :delta_x
        pmap[:Upsilon] = true
        dUpsilon[:, :] .= 0.
        dUpsilon[J+1, J+2] = 1.
    ##### elasticity params
    elseif param_sym == :a0
        pmap[:alpha0] = true
        pmap[:alpha0_til] = true
        dalpha0[1] = 1.
        dalpha0_til[1] = 1.
    elseif param_sym == :n0
        pmap[:eta0] = true
        deta0[1] = 1.
    else
        throw(DomainError(param_sym, "param_sym not recognized in semiC model"))
    end
    return nothing
end

# derivs wrt AR_hat
function deriv_AR_semiC!(nkphsolvr::NKPHSolver, j_idx::Int)::Nothing
    # compute model object derivs wrt AR_hat_{j}
    J, N_params = nkphsolvr.J, nkphsolvr.N_params
    # pull out siginv parameter
    siginv = nkphsolvr.aux_macro_params[:siginv]
    deriv_idx = _AR_idx_to_deriv_idx(nkphsolvr, j_idx)
    dUpsilon = @view nkphsolvr.dUpsilon[:, :, deriv_idx]
    dUpsilon[:, :] .= 0.
    dUpsilon[J+2, j_idx] = -siginv
    return nothing
end

# init
function initialize_nkphsolvr_semiC(
        a::Float64,
        df_short_rate_mat::DataFrame,
        df_params_constraints::DataFrame,
        df_params_fixed::DataFrame,
    )::NKPHSolver
    # model objects
    ei_vec = Vector{Float64}(df_short_rate_mat[:, :ei_vec])
    ed_vec = Vector{Float64}(df_short_rate_mat[:, :ed_vec])
    param_names = [Symbol(pname) for pname in df_params_constraints[:, :param]]
    J = length(ei_vec)
    N = J + 2
    N_params = size(df_params_constraints, 1)

    # fixed params
    pfixed_dict = OrderedDict{Symbol, Float64}()
    pfixed_keys = [Symbol(pname) for pname in df_params_fixed[:, :param]]
    pfixed_vals = Vector{Float64}(df_params_fixed[:, :val])
    for (idx, key) in enumerate(pfixed_keys)
        pfixed_dict[key] = pfixed_vals[idx]
    end

    # create model objects with fixed params
    Theta0 = zeros(J)
    Theta0_til = zeros(J)
    Theta1_til = zeros(J)
    Theta0[5] = pfixed_dict[:theta0_s]
    Theta0[6] = pfixed_dict[:theta0_l]
    Theta1 = zeros(J)
    Theta1[5] = pfixed_dict[:theta1_s]
    Theta1[6] = pfixed_dict[:theta1_l]
    Theta0_til[7] = pfixed_dict[:theta0_til]
    Theta1_til[7] = pfixed_dict[:theta1_til]

    # elasticities
    a1 = pfixed_dict[:alpha1]
    a1_til = pfixed_dict[:alpha1_til]
    n1 = pfixed_dict[:eta1_til]

    # placeholders
    sigma = zeros(J, J)
    a0 = 0.
    a0_til = 0.
    n0 = 0.
    # solver object
    nkphsolvr = NKPHSolver(
        map_params_semiC!, deriv_param_semiC!, deriv_AR_semiC!,
        param_names,
        N, J,
        sigma, a, a0, a1, a0_til, a1_til, n0, n1,
        Theta0, Theta1, Theta0_til, Theta1_til,
        ei_vec, ed_vec,
        N_params
    )
    # add aux macro variables
    nkphsolvr.aux_macro_params[:siginv] = pfixed_dict[:siginv]
    nkphsolvr.aux_macro_params[:rho_pi] = pfixed_dict[:rho_pi]
    return nkphsolvr
end



#############################################################################################
#############################################################################################
# semi-simplified "C" QE model
# no demand factor reaction to macro variables
# demand factor reaction to short rate
# Treasury and Risky demand factors identical parameterization
# alpha and alpha_til identical parameterization
# siginv and rho_pi fixed
# single QE factor

function map_params_semiC_QE!(nkphsolvr::NKPHSolver)::Nothing
    # update all model objects
    # pull out objects
    M, Upsilon, J, N_params = nkphsolvr.M, nkphsolvr.Upsilon, nkphsolvr.J, nkphsolvr.N_params
    sdyn = nkphsolvr.sdyn
    eta0, eta1 = nkphsolvr.eta0, nkphsolvr.eta1
    sigma, a = nkphsolvr.sigma, nkphsolvr.a
    alpha0, alpha0_til = nkphsolvr.alpha0, nkphsolvr.alpha0_til
    AR_hat = nkphsolvr.AR_hat
    aux_macro_params = nkphsolvr.aux_macro_params

    # pull out individual params from x_dict
    x_dict = nkphsolvr.x_dict

    # update model objects (inplace)
    # direct updates
    alpha0[1] = x_dict[:a0]
    alpha0_til[1] = x_dict[:a0]
    eta0[1] = x_dict[:n0]

    sigma[1, 1] = x_dict[:s_i]
    sigma[2, 2] = x_dict[:s_d]
    sigma[3, 3] = x_dict[:s_zpi]
    sigma[4, 4] = x_dict[:s_zx]
    sigma[5, 5] = sigma[6, 6] = sigma[7, 7] = x_dict[:s_b]

    # dynamics
    # Upsilon matrix
    # i_t
    Upsilon[1, 1] = x_dict[:k_i]
    Upsilon[1, J+1] = -x_dict[:k_i] * x_dict[:p_ipi]
    Upsilon[1, J+2] = -x_dict[:k_i] * x_dict[:p_ix]
    # d_t
    Upsilon[2, 2] = x_dict[:k_d]
    Upsilon[2, J+1] = -x_dict[:k_d] * x_dict[:p_dpi]
    Upsilon[2, J+2] = -x_dict[:k_d] * x_dict[:p_dx]
    # z_pit
    Upsilon[3, 3] = x_dict[:k_zpi]
    # z_xt
    Upsilon[4, 4] = x_dict[:k_zx]
    # b_st/b_lt/btil_t
    Upsilon[5, 5] = Upsilon[6, 6] = Upsilon[7, 7] = x_dict[:k_b]
    #Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:k_b] * x_dict[:p_bi]
    Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:p_bi]
    # pi_t
    Upsilon[J+1, 3] = 1.
    Upsilon[J+1, J+1] = -aux_macro_params[:rho_pi]
    Upsilon[J+1, J+2] = x_dict[:delta_x]
    # x_t
    Upsilon[J+2, 1:J] = -aux_macro_params[:siginv] .* AR_hat
    Upsilon[J+2, 4] += aux_macro_params[:siginv]

    # additional QE variables
    Upsilon[J, J] = nkphsolvr.aux_macro_params[:kappa_QE]
    sigma[J, J] = nkphsolvr.aux_macro_params[:sigma_QE]
    return nothing
end

# init
function initialize_nkphsolvr_semiC_QE(
        a::Float64,
        df_short_rate_mat::DataFrame,
        df_params_constraints::DataFrame,
        df_params_fixed::DataFrame,
    )::NKPHSolver
    # model objects
    ei_vec = Vector{Float64}(df_short_rate_mat[:, :ei_vec])
    ed_vec = Vector{Float64}(df_short_rate_mat[:, :ed_vec])
    param_names = [Symbol(pname) for pname in df_params_constraints[:, :param]]
    J = length(ei_vec)
    N = J + 2
    N_params = size(df_params_constraints, 1)

    # fixed params
    pfixed_dict = OrderedDict{Symbol, Float64}()
    pfixed_keys = [Symbol(pname) for pname in df_params_fixed[:, :param]]
    pfixed_vals = Vector{Float64}(df_params_fixed[:, :val])
    for (idx, key) in enumerate(pfixed_keys)
        pfixed_dict[key] = pfixed_vals[idx]
    end

    # create model objects with fixed params
    Theta0 = zeros(J)
    Theta0_til = zeros(J)
    Theta1_til = zeros(J)
    Theta0[5] = pfixed_dict[:theta0_s]
    Theta0[6] = pfixed_dict[:theta0_l]
    Theta1 = zeros(J)
    Theta1[5] = pfixed_dict[:theta1_s]
    Theta1[6] = pfixed_dict[:theta1_l]
    Theta0_til[7] = pfixed_dict[:theta0_til]
    Theta1_til[7] = pfixed_dict[:theta1_til]

    # additional QE factors
    Theta0[8] = pfixed_dict[:theta0_QE]
    Theta1[8] = pfixed_dict[:theta1_QE]
    Theta0_til[8] = pfixed_dict[:theta0_til_QE]
    Theta1_til[8] = pfixed_dict[:theta1_til_QE]

    # elasticities
    a1 = pfixed_dict[:alpha1]
    a1_til = pfixed_dict[:alpha1_til]
    n1 = pfixed_dict[:eta1_til]

    # placeholders
    sigma = zeros(J, J)
    a0 = 0.
    a0_til = 0.
    n0 = 0.
    # solver object
    nkphsolvr = NKPHSolver(
        map_params_semiC_QE!, deriv_param_semiC!, deriv_AR_semiC!,
        param_names,
        N, J,
        sigma, a, a0, a1, a0_til, a1_til, n0, n1,
        Theta0, Theta1, Theta0_til, Theta1_til,
        ei_vec, ed_vec,
        N_params
    )
    # add aux macro variables
    nkphsolvr.aux_macro_params[:siginv] = pfixed_dict[:siginv]
    nkphsolvr.aux_macro_params[:rho_pi] = pfixed_dict[:rho_pi]
    # additional QE variables
    nkphsolvr.aux_macro_params[:kappa_QE] = pfixed_dict[:kappa_QE]
    nkphsolvr.aux_macro_params[:sigma_QE] = pfixed_dict[:sigma_QE]
    return nkphsolvr
end





#############################################################################################
#############################################################################################
# semi-simplified "C" QE passive model
# no demand factor reaction to macro variables
# demand factor reaction to short rate
# Treasury and Risky demand factors identical parameterization
# alpha and alpha_til identical parameterization
# siginv and rho_pi fixed
# two factor QE

function map_params_semiC_QE_passive!(nkphsolvr::NKPHSolver)::Nothing
    # update all model objects
    # pull out objects
    M, Upsilon, J, N_params = nkphsolvr.M, nkphsolvr.Upsilon, nkphsolvr.J, nkphsolvr.N_params
    sdyn = nkphsolvr.sdyn
    eta0, eta1 = nkphsolvr.eta0, nkphsolvr.eta1
    sigma, a = nkphsolvr.sigma, nkphsolvr.a
    alpha0, alpha0_til = nkphsolvr.alpha0, nkphsolvr.alpha0_til
    AR_hat = nkphsolvr.AR_hat
    aux_macro_params = nkphsolvr.aux_macro_params

    # pull out individual params from x_dict
    x_dict = nkphsolvr.x_dict

    # update model objects (inplace)
    # direct updates
    alpha0[1] = x_dict[:a0]
    alpha0_til[1] = x_dict[:a0]
    eta0[1] = x_dict[:n0]

    sigma[1, 1] = x_dict[:s_i]
    sigma[2, 2] = x_dict[:s_d]
    sigma[3, 3] = x_dict[:s_zpi]
    sigma[4, 4] = x_dict[:s_zx]
    sigma[5, 5] = sigma[6, 6] = sigma[7, 7] = x_dict[:s_b]

    # dynamics
    # Upsilon matrix
    # i_t
    Upsilon[1, 1] = x_dict[:k_i]
    Upsilon[1, J+1] = -x_dict[:k_i] * x_dict[:p_ipi]
    Upsilon[1, J+2] = -x_dict[:k_i] * x_dict[:p_ix]
    # d_t
    Upsilon[2, 2] = x_dict[:k_d]
    Upsilon[2, J+1] = -x_dict[:k_d] * x_dict[:p_dpi]
    Upsilon[2, J+2] = -x_dict[:k_d] * x_dict[:p_dx]
    # z_pit
    Upsilon[3, 3] = x_dict[:k_zpi]
    # z_xt
    Upsilon[4, 4] = x_dict[:k_zx]
    # b_st/b_lt/btil_t
    # TODO TODO TODO
    Upsilon[5, 5] = Upsilon[6, 6] = Upsilon[7, 7] = x_dict[:k_b]
    #Upsilon[5, 5] = x_dict[:k_b]
    #Upsilon[6, 6] = x_dict[:k_b] + 1e-6
    #Upsilon[7, 7] = x_dict[:k_b] + 2e-6

    Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:p_bi]
    # pi_t
    Upsilon[J+1, 3] = 1.
    Upsilon[J+1, J+1] = -aux_macro_params[:rho_pi]
    Upsilon[J+1, J+2] = x_dict[:delta_x]
    # x_t
    Upsilon[J+2, 1:J] = -aux_macro_params[:siginv] .* AR_hat
    Upsilon[J+2, 4] += aux_macro_params[:siginv]

    # additional QE variables (active and passive)
    Upsilon[J-1, J-1] = nkphsolvr.aux_macro_params[:kappa_QE_p]
    Upsilon[J, J] = nkphsolvr.aux_macro_params[:kappa_QE]
    Upsilon[J-1, J] = -nkphsolvr.aux_macro_params[:gamma_QE_p]
    sigma[J-1, J-1] = nkphsolvr.aux_macro_params[:sigma_QE_p]
    sigma[J, J] = nkphsolvr.aux_macro_params[:sigma_QE]
    return nothing
end

# init
function initialize_nkphsolvr_semiC_QE_passive(
        a::Float64,
        df_short_rate_mat::DataFrame,
        df_params_constraints::DataFrame,
        df_params_fixed::DataFrame,
    )::NKPHSolver
    # model objects
    ei_vec = Vector{Float64}(df_short_rate_mat[:, :ei_vec])
    ed_vec = Vector{Float64}(df_short_rate_mat[:, :ed_vec])
    param_names = [Symbol(pname) for pname in df_params_constraints[:, :param]]
    J = length(ei_vec)
    N = J + 2
    N_params = size(df_params_constraints, 1)

    # fixed params
    pfixed_dict = OrderedDict{Symbol, Float64}()
    pfixed_keys = [Symbol(pname) for pname in df_params_fixed[:, :param]]
    pfixed_vals = Vector{Float64}(df_params_fixed[:, :val])
    for (idx, key) in enumerate(pfixed_keys)
        pfixed_dict[key] = pfixed_vals[idx]
    end

    # create model objects with fixed params
    Theta0 = zeros(J)
    Theta0_til = zeros(J)
    Theta1_til = zeros(J)
    Theta0[5] = pfixed_dict[:theta0_s]
    Theta0[6] = pfixed_dict[:theta0_l]
    Theta1 = zeros(J)
    Theta1[5] = pfixed_dict[:theta1_s]
    Theta1[6] = pfixed_dict[:theta1_l]
    Theta0_til[7] = pfixed_dict[:theta0_til]
    Theta1_til[7] = pfixed_dict[:theta1_til]

    # additional QE factors
    Theta0[8] = pfixed_dict[:theta0_QE_p]
    Theta0[9] = pfixed_dict[:theta0_QE]
    Theta0_til[8] = pfixed_dict[:theta0_til_QE_p]
    Theta0_til[9] = pfixed_dict[:theta0_til_QE]

    Theta1[8] = pfixed_dict[:theta1_QE]
    Theta1[9] = pfixed_dict[:theta1_QE]
    Theta1_til[8] = pfixed_dict[:theta1_til_QE]
    Theta1_til[9] = pfixed_dict[:theta1_til_QE]

    # elasticities
    a1 = pfixed_dict[:alpha1]
    a1_til = pfixed_dict[:alpha1_til]
    n1 = pfixed_dict[:eta1_til]

    # placeholders
    sigma = zeros(J, J)
    a0 = 0.
    a0_til = 0.
    n0 = 0.
    # solver object
    nkphsolvr = NKPHSolver(
        map_params_semiC_QE_passive!, deriv_param_semiC!, deriv_AR_semiC!,
        param_names,
        N, J,
        sigma, a, a0, a1, a0_til, a1_til, n0, n1,
        Theta0, Theta1, Theta0_til, Theta1_til,
        ei_vec, ed_vec,
        N_params
    )
    # add aux macro variables
    nkphsolvr.aux_macro_params[:siginv] = pfixed_dict[:siginv]
    nkphsolvr.aux_macro_params[:rho_pi] = pfixed_dict[:rho_pi]
    # additional QE variables
    nkphsolvr.aux_macro_params[:kappa_QE] = pfixed_dict[:kappa_QE]
    nkphsolvr.aux_macro_params[:kappa_QE_p] = pfixed_dict[:kappa_QE_p]
    nkphsolvr.aux_macro_params[:gamma_QE_p] = pfixed_dict[:gamma_QE_p]
    nkphsolvr.aux_macro_params[:sigma_QE] = pfixed_dict[:sigma_QE]
    nkphsolvr.aux_macro_params[:sigma_QE_p] = pfixed_dict[:sigma_QE_p]
    return nkphsolvr
end







#############################################################################################
#############################################################################################
# semi-simplified "C" MP model
# no demand factor reaction to macro variables
# demand factor reaction to short rate
# Treasury and Risky demand factors identical parameterization
# alpha and alpha_til identical parameterization
# siginv and rho_pi fixed

function map_params_semiC_MP!(nkphsolvr::NKPHSolver)::Nothing
    # update all model objects
    # pull out objects
    M, Upsilon, J, N_params = nkphsolvr.M, nkphsolvr.Upsilon, nkphsolvr.J, nkphsolvr.N_params
    sdyn = nkphsolvr.sdyn
    eta0, eta1 = nkphsolvr.eta0, nkphsolvr.eta1
    sigma, a = nkphsolvr.sigma, nkphsolvr.a
    alpha0, alpha0_til = nkphsolvr.alpha0, nkphsolvr.alpha0_til
    AR_hat = nkphsolvr.AR_hat
    aux_macro_params = nkphsolvr.aux_macro_params

    # pull out individual params from x_dict
    x_dict = nkphsolvr.x_dict

    # update model objects (inplace)
    # direct updates
    alpha0[1] = x_dict[:a0]
    alpha0_til[1] = x_dict[:a0]
    eta0[1] = x_dict[:n0]

    sigma[1, 1] = x_dict[:s_i]
    sigma[2, 2] = x_dict[:s_d]
    sigma[3, 3] = x_dict[:s_zpi]
    sigma[4, 4] = x_dict[:s_zx]
    sigma[5, 5] = sigma[6, 6] = sigma[7, 7] = x_dict[:s_b]

    # dynamics
    # Upsilon matrix
    # i_t
    Upsilon[1, 1] = x_dict[:k_i]
    Upsilon[1, J+1] = -x_dict[:k_i] * x_dict[:p_ipi]
    Upsilon[1, J+2] = -x_dict[:k_i] * x_dict[:p_ix]
    # d_t
    Upsilon[2, 2] = x_dict[:k_d]
    Upsilon[2, J+1] = -x_dict[:k_d] * x_dict[:p_dpi]
    Upsilon[2, J+2] = -x_dict[:k_d] * x_dict[:p_dx]
    # z_pit
    Upsilon[3, 3] = x_dict[:k_zpi]
    # z_xt
    Upsilon[4, 4] = x_dict[:k_zx]
    # b_st/b_lt/btil_t
    Upsilon[5, 5] = Upsilon[6, 6] = Upsilon[7, 7] = x_dict[:k_b]
    #Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:k_b] * x_dict[:p_bi]
    Upsilon[5, 1] = Upsilon[6, 1] = Upsilon[7, 1] = -x_dict[:p_bi]
    # pi_t
    Upsilon[J+1, 3] = 1.
    Upsilon[J+1, J+1] = -aux_macro_params[:rho_pi]
    Upsilon[J+1, J+2] = x_dict[:delta_x]
    # x_t
    Upsilon[J+2, 1:J] = -aux_macro_params[:siginv] .* AR_hat
    Upsilon[J+2, 4] += aux_macro_params[:siginv]

    # additional MP variables
    Upsilon[J, J] = nkphsolvr.aux_macro_params[:kappa_MP]
    Upsilon[5, J] = Upsilon[6, J] = Upsilon[7, J] = -x_dict[:p_bi]
    return nothing
end

# init
function initialize_nkphsolvr_semiC_MP(
        a::Float64,
        df_short_rate_mat::DataFrame,
        df_params_constraints::DataFrame,
        df_params_fixed::DataFrame,
    )::NKPHSolver
    # model objects
    ei_vec = Vector{Float64}(df_short_rate_mat[:, :ei_vec])
    ed_vec = Vector{Float64}(df_short_rate_mat[:, :ed_vec])
    param_names = [Symbol(pname) for pname in df_params_constraints[:, :param]]
    J = length(ei_vec)
    N = J + 2
    N_params = size(df_params_constraints, 1)

    # fixed params
    pfixed_dict = OrderedDict{Symbol, Float64}()
    pfixed_keys = [Symbol(pname) for pname in df_params_fixed[:, :param]]
    pfixed_vals = Vector{Float64}(df_params_fixed[:, :val])
    for (idx, key) in enumerate(pfixed_keys)
        pfixed_dict[key] = pfixed_vals[idx]
    end

    # create model objects with fixed params
    Theta0 = zeros(J)
    Theta0_til = zeros(J)
    Theta1_til = zeros(J)
    Theta0[5] = pfixed_dict[:theta0_s]
    Theta0[6] = pfixed_dict[:theta0_l]
    Theta1 = zeros(J)
    Theta1[5] = pfixed_dict[:theta1_s]
    Theta1[6] = pfixed_dict[:theta1_l]
    Theta0_til[7] = pfixed_dict[:theta0_til]
    Theta1_til[7] = pfixed_dict[:theta1_til]

    # elasticities
    a1 = pfixed_dict[:alpha1]
    a1_til = pfixed_dict[:alpha1_til]
    n1 = pfixed_dict[:eta1_til]

    # placeholders
    sigma = zeros(J, J)
    a0 = 0.
    a0_til = 0.
    n0 = 0.
    # solver object
    nkphsolvr = NKPHSolver(
        map_params_semiC_MP!, deriv_param_semiC!, deriv_AR_semiC!,
        param_names,
        N, J,
        sigma, a, a0, a1, a0_til, a1_til, n0, n1,
        Theta0, Theta1, Theta0_til, Theta1_til,
        ei_vec, ed_vec,
        N_params
    )
    # add aux macro variables
    nkphsolvr.aux_macro_params[:siginv] = pfixed_dict[:siginv]
    nkphsolvr.aux_macro_params[:rho_pi] = pfixed_dict[:rho_pi]
    # additional MP variables
    nkphsolvr.aux_macro_params[:kappa_MP] = pfixed_dict[:kappa_MP]
    return nkphsolvr
end
